The Origin of the Fragile-to-Strong Crossover in Liquid Silica as Expressed by its 

Potential Energy Landscape 
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The origin of the fragile-to-strong crossover in liquid silica is characterized in terms of properties 
of the potential energy landscape (PEL). Using the standard BKS model of silica we observe a low- 
energy cutoff of the PEL. It is shown that this feature of the PEL is responsible for the occurrence of 
the fragile-to-strong crossover and may also explain the avoidance of the Kauzmann paradox. The 
number of defects, i.e. deviations from the ideal tetrahedral structure, vanishes for configurations 
with energies close to this cutoff. This suggests a structural reason for this cutoff. 

PACS numbers: 64.70.Pf,65.40.Gr,66.20.+d 



Understanding the properties of liquid silica is of ut- 
most importance due to its enormous technological and 
scientific relevance. In the Angell [l| classification scheme 
it is a very strong glass, i.e. displaying an Arrhcnius tem- 
perature dependence of its transport observables like vis- 
cosity or difFusivity. Furthermore one may expect that 
as a typical network-forming system it shares properties 
with other network-formers like water. Indeed, recently it 
has been pointed out that anomalies like negative thermal 
expansion, isobaric heat capacity minima and isothermal 
compressibility minima occur in water as well as in sim- 
ulated liquid silicaHHQ- 

The experimental oxygen diffusion constant D(T) be- 
haves like D exp(-E a /k B T) with D = 2.6 cm 2 /s 0. 
For a microscopic hopping time of tq = 10~ 13 s and a 
typical length scale of &—2 A one would expect Dq w 
a 2 /(6ro) ~ 10~ 3 cm 2 /s. This discrepancy of some orders 
of magnitude could be rationalized if for temperatures 
above the region, which is accessible in most experiments, 
the diffusivity D(T) would show deviations from the Ar- 
rhcnius behavior. Indeed, from experiments 6] as well 
as from simulations Q this so-called "fragile-to-strong" 
crossover (FSC) has been found for silica in a temper- 
ature range around Tpsc ~ 3500 K. A similar scenario 
has been discussed for water Q and BeF 2 0- It has 
been speculated that the occurrence of a FSC is directly 
related to the presence of polyamorphism 0, ^| . 

In recent years the analysis of the potential energy 
landscape (PEL) became an indispensable tool to eluci- 
date the thermodynamics and dynamics of many different 
supercooled liquids Ell IS 13, EE IS llE 113, 
l2l[ I22I |23| . The dynamics of the total system is repre- 
sented as the dynamics of a point in the high-dimensional 
configurational space. In particular, the thermodynam- 
ics can be expressed in terms of the distribution of the 
inherent structures (minima of the PEL) [l4l |24| . It 
was also possible to express the long-range dynamics, i.e. 
the diffusivity, in terms of PEL-characteristics |20l |25j . 
Thus the phenomenon of the FSC should be reflected 
by properties of the PEL. Actually, for liquid silica it 



has been shown that around Tp$c the temperature de- 
pendence of the average inherent structure energy and, 
correspondingly, the configurational entropy S C (T) shows 
an inflection point 1 1 11 ]- This has direct consequences 
for the Kauzmann paradoxon. Extrapolating the S C (T)- 
dependence for T > Tpsc to lower temperatures one 
would obtain S c (Tk) — for a finite temperature Tk, 
the Kauzmann temperature. The inflection point around 
T w Tpsc suggests that S C (T) > for all temperatures. 
Thus it is speculated that "the FSC and the Kauzmann 
paradox are fundamentally interrelated phenomena" . 
In this work we will show that this is indeed the case and 
suggest the underlying reason for both properties. 

For the modelling of silica we use the standard BKS 
potential j2(J. As compared to the standard choice we 
employ a somewhat larger cutoff-radius of 8.5 A for the 
short-range interaction to avoid energy drifts. Simula- 
tions are performed in the NVE ensemble with periodic 
boundary conditions, using a density of p = 2.30 g/cm 3 . 
The lowest simulated temperature was 2800 K. In re- 
cent years the concept of metabasins (MB) has been in- 
troduced as an appropriate way of coarse-graining the 
configuration space [23, 12H |29| . On a qualitative level, 
adjacent inherent structures between which the system 
performs back- and forth jumps are merged together to 
one MB (see [2^| for a precise definition). The energy e 
of a MB is defined as the energy of the inherent structure 
with the lowest energy. Whereas for the thermodynamics 
the distinction between inherent structures and MBs is 
not relevant, the relation to the dynamics is much simpler 
in terms of MBs. For a binary Lennard-Jones system it 
could be shown that the average residence time in a MB 
but not that in an inherent structure is directly propor- 
tional to the inverse diffusion constant [25|. 

To characterize the PEL it is crucial to analyse rather 
small systems with periodic boundary conditions [22ll3fj| . 
Otherwise, interesting aspects may be simply averaged 
out, as shown below. For this purpose one may look 
for the minimum system size without any significant fi- 
nite size effects. For a binary Lennard-Jones system it 
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could be shown that starting from circa 65 particles no 
finite-size effects are present for the ran ge o f tempera- 
tures accessible to computer simulations 31]. For BKS 
silica it is known that the diffusion constant becomes con- 
stant only for system sizes of more than 1000 particles 
[13 • We have calculated the incoherent scattering func- 
tion S(q ma x,t) for system sizes ranging from 60 to 1000 
at T T c . It turns out that they all have identical non- 
exponential relaxation behavior in the a-regime except 
for simple scaling (e.g. r a {N = 99)/r Q (iV = 1000) w 3); 
see inset of Fig.l. This is a first hint that the transport 
behavior is identical for all system sizes in this range. In 
what follows we report simulations for N = 99 (denoted 
as BKS99). 

The temperature dependence of the oxygen diffusion 
constant D(T) of BKS99 is shown in Fig.la. It dis- 
plays the same FSC as the macroscopic system with Ar- 
rhenius behavior below 3500 K and a very similar low- 
temperature activation energy (4.9 eV). In Fig. lb we 
show the average energy of MBs (e(T)) and inherent 
structures (e/g(T)). Furthermore we include the inherent 
structures energies from Ref. [Tlj where a much larger 
BKS silica system (N = 1332) with slightly different cut- 
off conditions and nearly the same density (p = 2.31 
gem -3 ) has been simulated. The temperature depen- 
dence is nearly identical for all three curves and has 
a change in slope around Tpsc- We note in passing 
that the temperature dependence of the dynamic het- 
erogeneities in real space does not show any anomalies 
around Tpsc Wk- 

So far we have reproduced the temperature dependence 
of the dynamics (FSC) and thermodynamics (inflection 
in (e(T)), related to the Kauzmann paradox) for BKS99. 
Going beyond previous work the main goal of this work 
is to analyse the relevant observables in terms of energy. 
From this we will uniquely identify a common origin of 
the FSC and the inflection in (e(T)): just around Tpsc 
the system approaches the low-energy end of the PEL. 

One of the most fundamental energy-dependent quan- 
tities is the density of MBs G(e). In a first step 
we record the distribution of MBs p(e, T) at differ- 
ent temperatures. In case that anharmonicities are 
small or do not strongly depend on temperature one 
has p(e,T) cx G(e) exp(-/3e) with 3 = l/(k B T) or, 
equivalents, G{e) cx p(e, T) exp(/3e) [3II3- Thus via 
reweighting it is possible to extract G(e) from p(e, T) 
apart from a proportionality constant. Actually, for 
T > 5000 K this reweighting breaks down because of 
anharmonicities. We have estimated the proportionality 
constant from the condition that for T = 5000 K the con- 
figurational entropy is identical to the value of 6.8 J/(mol 
K), reported in [111 ]. The resulting energy-distribution 
G(e) is shown in Fig. 2. It is possible to fit G(e) with 
a gaussian G(e) 
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to the binary Lennard-Jones system the distribution is 
gaussian within statistical errors. 

Naturally, during our simulations we find a lowest MB 
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Figure 1: (a) Temperature dependence of the oxygen diffu- 
sion constant for BKS99. In the inset the incoherent scatter- 
ing function S(q ma x = 1.7A _1 ,t) is shown at T = 3300 K 
and 4000 K for system sizes between 60 and 1000 and after 
scaling along both axes, (b) Temperature dependence of the 
average energy of MBs (e(T)), the average energy of inherent 
structures (eis(T)). Included are data from Ref.[ll|, shifted 
by 0.06 eV. The broken line describes the gaussian approxi- 
mation with parameters obtained from Fig. 2. 



energy, which we denote as e c (e c = —1910.7 eV). In 
principle there are two possible scenarios for this cut- 
off: (1) The density of states G(e) is also Gaussian for 
e < e c . These states, however, are not found due to finite 
simulation times. The resulting distribution is denoted 
G\(e). (2) The number of states below e c is much smaller 
than predicted by the Gaussian extrapolation such that 
e c can be basically regarded as a low-energy cutoff of the 
energy landscape (denoted G cu toff(e))- To falsify sce- 
nario (1) we assume the validity of G\{e). In a first step 
one can compare the average ener gy with the gaussian 
prediction (e(T)) = e - a 2 /T [llll|. It is included in 
Fig. lb. One can clearly see the deviations for T < 4000 
K, indicating that the description by a gaussian breaks 
down at low energies. Actually, the deviations at high 
temperatures are due to anharmonicities. An even closer 
analysis can be performed by comparing the whole dis- 
tribution p(e,T), obtained from our simulations, with 
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Figure 2: The energy distribution G(e), obtained from 
reweighting of p(e, T) for temperatures between 2800 K and 
5000 K. The resulting gaussian-type fits Gi(e) and G2(e) are 
included (see text). 



Pi(e,T) oc G\{e) exp(— /3e). This comparison is shown 
in Fig. 3. The dramatic difference between pi(e,T=2800 
K) and p(e, T=2800 K) again shows that the gaussian 
extrapolation is not valid. Actually, this difference can 
be also expressed by stating that more than 100 MBs 
with energies below e c should have been found during 
our simulation time of 120 ns at T = 2800 K in case 
of scenario (1). This quantification, however, is rather 
involved. Since the discrepancy is already obvious from 
Fig. 3 we defer it to a later publication. In any event, we 
have to conclude that the appearance of a cutoff is not an 
artifact but reflects a major depletion of states below e c 
as compared to the gaussian distribution. Actually, using 
the classification scheme recently introduced by Ruocco 
et al. this result clearly shows that BKS silica is a "B 
strong liquid" 34] . 

Before validating scenario (2) we slightly smear out 
the cutoff by introducing the distribution C?2(e) oc 
Gcutoff(e) ■ tanh((e — e c ) 2 /2), defined for e > e c (see 
Figj^J. As shown in Fig. 3 the estimation p2(e,T) oc 
G*2(e) exp(— j3e) reproduces p(e, T) for all temperatures. 

This analysis may help to clarify why it is so important 
to use rather small systems. Based on the assumption 
that a system of n ■ 99 (ri = 2, 3) particles is composed of 
n independent but identical subsets of 99 particles it pos- 
sible to estimate p(e/n,T = 2800 K); see inset of Fig. 3. 
Whereas for 99 particles the cutoff can be directly seen 
from the strongly asymmetric shape of p(e, T), this ef- 
fect is smeared out for 198 or 297 particles and in the 
limit of large n, p(e/n, T) would finally approach a delta- 
function. Thus the presence of the cutoff in any subsys- 
tem is no longer visible for larger systems. This line of 
thought explicitly rationalizes the use of small systems. 

Whereas for high temperatures the thermodynamics 
and the dynamics is not influenced by the cutoff at e c , 
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Figure 3: The simulated distribution p(e, T) for T=2800 K, 
T=3500 K, and T=4000 K. Furthermore we include Pl (e) 
and p2(e,T) (for the same three temperatures). Note that 
Pi(e,T = 4000 K) » p 2 (e,T = 4000 K). In the inset the 
estimation of p2(e/n, T = 28007?) for different systems sizes 
n ■ 99 (n=l,2,3) is displayed. 



a significant impact is expected at low temperatures. As 
seen from Fig. 3 for temperatures below 4000 K the cutoff- 
population p(e w e Cl T) starts to become relevant, re- 
sulting in an increasing difference between p\{e,T) and 
p(e,T). This is exactly the temperature range of the 
FSC. We suggest that this agreement is not coincidental. 
As shown for the more fragile binary Lennard- Jones sys- 
tem the increase of the apparent activation energy with 
decreasing temperature can be related to the fact that 
the system explores lower regions of the PEL at low tem- 
peratures for which the local effective saddles to escape 
the MB are correspondingly higher 22]. Not surpris- 
ingly, the same phenomenon is seen for BKS99 (data not 
shown). Thus the growing influence of the bottom of the 
PEL with decreasing temperature stops the increase of 
the apparent activation energy. 

What is the microscopic origin of the low-energy cut- 
off? One may expect that the energy of a structure is 
correlated with the number of defects, i.e. silicon atoms 
with three or five oxygens in the nearest neighbor shell 
or oxygens with one or three silicon atoms. This is sug- 
gested by the significant decrease of the number of defects 
with decreasing temperature, observed in previous simu- 
lations 0, El3 • I n a nrs t step, we compare the occurrence 
of defects with the results for BKS8000 ■ It turns out 
that the number of different defects of the equilibrium 
configurations at T = 4000 K on average agree within a 
factor of 1.15 with the numbers, reported in 0. Thus 
the properties of defects do not show major finite size 
effects. In a second step we determine the fraction of 
these defects (nearest neighbor-shell: rsio < 2.28 A) for 
the minimized structures in dependence of energy; see 
Fig. 4. One can clearly see that the fraction of defects 
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Figure 4: The fraction of silicon (oxygen) atoms with 3 or 
5 oxygen (1 or 3 silicon) atoms in the nearest neighbor shell 
for different temperatures (right axis; silicon: stars; oxygen: 
squares) and the width of the partial Si-O structure factor for 
the nearest-neighbor shell (left axis). 



varies roughly linearly with energy and approaches zero 
for configurations around e ~ e c . These configurations, 
however, are still amorphous. This is checked by cal- 
culating the width of the first nearest-neighbor peak of 
the Si-0 partial structure factor for the minimized struc- 
tures. It is a measure for the degree of disorder in the 
nearest-neighbor shell. This value only weakly depends 
on energy; see Fig. 4. We may conclude that the reduc- 
tion of the number of defects is an efficient way to obtain 
configurations with lower energy. After reaching con- 
figurations with no defects at all there emerges a lower 
energy limit which cannot be crossed without crystalli- 



sation. This suggests a structural reason for this cutoff. 
In particular this excludes the scenario that the cutoff is 
a trivial consequence of the fact that for a finite number 
of MBs (or inherent structures) at some temperature the 
lowest energy configuration becomes thermodynamically 
relevant. As seen in Fig. 2 we indeed estimate an expo- 
nentially large number of different states around e ~ e c . 

This scenario is very different as compared to the frag- 
ile binary Lennard- Jones system where pi(e, T) ps p(e, T) 
for all simulated temperatures [2{| , excluding a PEL cut- 
off in the accessible energy range. Without network con- 
straints the system can optimize the disordered structure 
more efficiently to find configurations with even lower en- 
ergies (see also [35)1. 

Is it possible that the FSC for macroscopic BKS is of 
different nature than for BKS99? One might think that 
for larger systems it is possible to form more low-energy 
states due to the smaller number of constraints related 
to the periodic boundary conditions. Then G(e) might 
extend the gaussian-like behaviour also to lower energies 
and the temperature dependence of (e(T)) should ap- 
proach the gaussian limit. Comparison of BKS99 with 
BKS1332 shows that this effect, if at all, is very small 
(Fig.l). Therefore macroscopic BKS not only displays 
the same phenomenology as BKS99 but very likely also 
follows similar underlying mechanisms. 

In summary, for the prototype network former silica we 
have related the fragile-to-strong crossover as well as the 
inflection of (e(T)) to the presence of a lower cutoff in the 
PEL. Since this cutoff is related to structural constraints 
for disordered structures of a network former, it may be 
of general nature. 
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